Coexisting ordinary elasticity and superfluidity in a model of defect-free supersolid 
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We present the mechanics of a model of supersolid in the frame of the Gross- 
Pitaevskii equation at T = OK that do not require defects nor vacancies. A set 
of coupled nonlinear partial differential equations plus boundary conditions is 
derived. The mechanical equilibrium is studied under external constrains as 
steady rotation or external stress. Our model displays a paradoxical behavior: 
the existence of a non classical rotational inertia fraction in the limit of small 
rotation speed and no superflow under small (but finite) stress nor external 
force. The only matter flow for finite stress is due to plasticity. 
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The recent surge of interest in supersolidsi makes it important to reach a clearer understanding of the mechanical 
properties of such materials. In particular why a supersolid behaviour is observed in a rotating experiment, whilst, 
as in a ordinary solid no constant mass flux is driven by pressure gradient? (see Refs^) In 3 - two of us (YP and SR) 
proposed a fully explicit model of supersolid where many properties can be discussed in details. We thought timely to 
reconsider this model, in particular with respect to its properties of elasticity coupled to its ability to carry some sort 
of superflow in the absence of any defect. Although supersolidity is often related to the presence of defects, vacancies 
and so forth our model introduces an important difference between ordinary (classical) crystals and supcrsolids: in 
perfect classical crystals there is an integer number (or a simple fraction) of atoms per unit cell. Therefore the number 
density and the lattice parameters are not independent. On the contrary, in our model of supersolid, there is no such 
relation. The lattice parameters and the average density can be changed independently. 

This model is based upon the original Gross-Pitaevksii equation (G-P later with an integral term with a kernel 
that can be seen as a two-body potential in a first Born approximation 1 . This model yields the exact spectrum found 
long ago by Bogoliubo^, namely a relation between the energy of the elementary excitations and their momentum de- 
pending on the two-body potential. In this framework the roton minimum is a precursor of crystallisation. Something 
predicted ir£ where the possibility of a linear instability was only considered, although the transition is subcritical- 
first order—. By increasing the density, crystallization happens through a first order phase transition (see figure [5]). 
As shown in 3 - the crystal phase shows a periodic modulation of density in space together with some superfluid-like 
behaviour under rotation. 

The aim of the present letter is to show that, besides this behaviour, the system has also solid-like behaviour, at 
least under small stress. At larger stress, it flows plastically, the plasticity being facilitated by the eventual presence 
of defects. We derive the equation of motion for the average density n, the phase $ and the displacement u in the 
solid. A new propagating mode appears besides the usual longitudinal and transverse phonons in regular crystals. 
This mode is partly a modulation of the coherent quantum phase, like the phonons in superfluids at zero temperature. 
We discuss at the end the boundary conditions and how to handle steady rotation and pressure driven flow in this 
model. 

Our starting point is the original G-P equation^ for the complex valued wavefunction i/j(r,t) common to all bosonic 
particles of mass m: 



m d 4 = -^ 

dt 2m 



V 2 V + V y dr'[/(r' ~r)\^(r',t)\ 2 , (1) 



where [/(•) is a two body potential depending on the distance. For the numerics we choose a potential U(\r\ 
Uo0(a — \r\), with 6(.) Heaviside function. 



1 This is an exact approximation for the full quantum problem in the van der Waals limit of long range low amplitude potential U (r) = 
7 3 W(7r) whenever 7^0. 
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The ground state is a solution of equation {T]) of the form i/>o(r)e -1 » . It is a crystal when ipoir) is a periodic 
function such that tpo(r + q a a + qbb + q c c) = Vo(t) f° r la, b, c arbitrary integers, a,b,c being vectors defining the 
elementary lattice cell. This solution is the ground state in the sense that, given an average number density n = 
^ J dr|^o| 2 , ^ being the total volume, the lattice parameters and the function ipo(r) make E$ the smallest possible. 
The local density n(r,t), the displacement field u(r,t) of the crystal lattice and the slowly varying phase $(r,t) of 

The Lagrangian density for the G-P equation {TJ) reads in polar coordinates, ip = -yfpe^: 

C = -hp^ - ^ (p(V0) 2 + ^(Vp) 2 ) - \ P {r) J dr'U(r' r)p{r'). 

The ground state is given by the solution of the nonlinear intcgro-differential equation for p derived by variation of 
the action taking the phase field <p uniform in space: <jy = —pt/H, p constant, that is a solution of: 

-M+^f^-^W<^-r),M = 0. (2) 



4m \ 2p 2 

This ground state solution, if periodic in space, as we shall assume in full agreement with our numerical results, 
depends on the dimensionless parameter A = Uo^^na 3 onlyi Although in Ref4 we discussed a ground state as a 
modulation close to an uniform density near the transition, that is for a finite roton gap, we have observed numerically 
that a crystal ground state exists in a wide range of densities. In the limit A ^> 1 the lattice tends to an array of 
sharp density pulses distant of a, the width of the pulse decreasing like A -1 / 2 . 

Let po(r\n) be a ground state solution, then po(r — u\n) is also a ground state solution with the same p, for a 
constant displacement field u. The general perturbations around the ground state allow that $, u and n become 
slow varying fields on space and time. As in Ref£ we follow the general method called homogenization^. This splits 
cleanly the long-wave behaviour of the various parameters and the short range periodic dependence upon the lattice 
parameters. Let write the Ansatz for density and phase: 

p = p (r - u,n(r,t)) + p(r - u, n, t) + . . . 

<t> = <f>(r,t) + 4>(r-u,n,t) + ... (3) 

where $, u and n are slow varying fields and 4> and p are small and fast varying periodic functions. Introducing this 
Ansatz into the Lagrangian @ one gets an effective Lagrangian made of four kind of terms: 

i) As n changes continuously the periodic solution of the intcgrodiffcrential equation say po{r), can be considered 
as a regular function of the Lagrange multiplier p imposing the average density n. Therefore by integrating over an 
unit cell (V) of the lattice the Lagrangian from which © is derived one obtains an averaged Lagrangian that depends 
on n only that writes 

- C n = £(n) = ^fJ v drpo(r) J U(r' - r)p (r')dr'. (4) 

This yields the simplest case of homogenization. 

ii) Similarly, terms mixing the slow varying phase field <&(/•, t) and po{f) can be averaged directly leading to 

£* = -n (hdt$ + |^r(V$) 2 ) where n = y J v p (r) dr. 

Next contributions need to solve the Eulcr-Lagrangc conditions for the fast variables <j> and p. We shall sketch the 
effective Lagrangian for the phase <f>, leaving for the reader a similar calculation for the deformation part. 

in) The <f> dependence term of this Lagrangian can be re- written of the form: = 

~^7i I ( 2 /°o^- ■ + po(V^) 2 ^ dr, where A = (V$ — (V<6 ■ V)« — j£d t u) (considered as a constant in the 
unit cell). The Euler-Lagrange condition for reads A ■ Vpo + V • (^po"V(pj = 0. This Poisson-like equation is to 

be solved within the unit cell of the lattice, for a function <j> that is periodic with the same period as pq. The result 
(that can be expressed as the minimum of a certain Rayleigh-Ritz functional) is linear in A and can be written 
as cf> = KiAi where K(r) is a vector-valued function of r that is periodic and satisfies Vipo + V • (po~VKi) = 0. 
Putting the result into the Lagrangian one obtains the relevant contribution for the slowly varying part of the 

phase: = ^QijAiAj with the positive defined matrix Qij = y J v po(r)^Ki ■ VKjdr. The Lagrange function 
depends on the slow variables only. We shall restrict ourselves below to crystal structures sufficiently symmetric 
to make diagonal Qij = g(n)Sij. The quantity g(n) is zero if the crystal modulation is absent and would be 
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very small for Bosc-Einstcin condensate with a non local interaction term. g(n) — > n when all the mass is strongly 
localized in the center of the cell site with a small overlap in between the different sites. This is presumably the 
situation of almost all materials in their solid state at low temperature. A large Young modulus is likely a measure 
of the small overlap of the wave functions from one site to the next, making 4 He exceptional at this respect. In other 
words when g(n) — > n the supersolid behaves as a ordinary solid state. 

iv) The same method of homogenization works for the long- wave perturbations of gradients of the displacement 
u and yields a contribution to the Lagrangian that reads C u = — |-^y/ci^ 1 ^ t - The coefficients Xijki are given by 
integrals over the unit cell of various functions defined explicitely. This is the familiar elastic energy of a Hookean 
solid. 

To summarize, the effective Lagrangian reads: 



eff 



-hn 



2m 



n(V$) - g(n) V<f> 



m Dm 



£{n) 



1 



dui duj~ 
dxj dxi 



(5) 



where ^ = ^ + • V«. This expression is remarkable because it is fully explicit for a given ground state of 

the G-P model. As one can check this Lagrangian is Galilean invariant. 

We conjecture that, because this Lagrangian satisfies the symmetries imposed by the underlying physics and because 
it includes a priori all terms with the right order of magnitude with respect to the derivatives, the general Lagrangian 
of any supersolid at zero temperature has the same structure. In a recent paper, Son£ derives a Galilean invariant 
Lagrangian such that (|5|) is a sub-class but with well defined coefficients like g(n), £(n) and \ju depending on the 
details of the crystal structure. 

The dynamical equations are derived by variation of the action J Cd 3 r dt which is seen as a functional of n, <!> and 
u. The variation with respect to n, u and $ yields (writing g'(n) = dg/dn, etc.): 



dt 



dt + 2m 
Dm 



g(n)( 



Dt 



(V$) 

m dxi ' 



V . ell \ , 1 V du i du k 

) +£{n)+ 2 X ^dx-d^ 
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dxj 



ijkl 



du k 
dxi 



dxi- 



Q( 



Dt 



S9# 3$ 

m dxi dxk 



% + -V ■ (nV$) - -d k [ g(n)(S lk - d k u t ) ( d& 
at m m 



m Duj 
~K ~Dt 



(6) 

(7) 

(8) 



The latter equation reduces to the familiar equation of mass conservation for potential flows whenever g(n) = 0, 
namely in the absence of modulation of the ground state. Although our equations of motion (|6I7I8|) and the one of 
Andrccv Lifshitz are almost identical in the zero temperature limit (see eqns. (16) of Ref^), our model has significant 
differences with their. Our solid cannot be considered as the normal component of a two "fluids" system, because it 
is on the same footing (phase coherent) as the superfluid part at T = OK. Therefore, at small finite temperature, 
our model has a normal component that is a fluid of vanishing density at T = OK, besides its coherent superfluid 
and solid part and should change the superfluid density fraction. Following Landau's ideas, this normal fluid is a gas 
of quasi-particles with the mixed spectrum able to carry momentum whilst the coherent part (superfluid plus solid) 
stays at rest. 

The Euler-Lagrange conditions impose also the boundary conditions for the equations of motion: 



nd k $ - g(Sik - d k ui) (d^ - 



efc = nVkSk- 



where Vk is the local speed of the solid wall of the container and Bk is normal to it. The displacement moves with the 
wall: ^ = V. 

Let us look at small perturbations around a nondeformed (u = 0) and steady (V$ = 0) state of average density n. 
The linearized version of (JGI7I8D shows that the shear waves are decoupled from the compression and phase (Bogoliubov- 
like) waves. The dispersion relation for the coupled compression and phase waves leads to a simple algebraic equation. 
In the limit g(n) — > the crystal structure disappears and the phase mode propagates at the usual speed of sound 

that is whenever the supersolid behaves as a 

+ c 2 and 



K 



found by Bogoliubov c = y/£"(n)/(mn). In the limit g(n) 

regular solid state, the two propagation speeds are [ck is the longitudinal elastic wave speech) v\ = y< 
?j 2 = y 1 c\c 2 /{c 2 K + c 2 )\J\ — g(n)/n meaning that the phase mode disappears at the transition supersolid-solid. 

As suggested by Leggetljii an Andronikashvili kind of experiment could manifest a non classical rotational inertia 
(NCRI). Indeed let us supose that the wall of the container of volume il rotate with an uniform angular speed u>. 
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Then for low angular speed the crystal moves rigidly with the container ii = u> x r without any elastic deformation. 
The densities n and g(n) being constant in space, equation ((SJ) simplifies into 

V 2 $ = in with V$ • e = (m/h)(u x r) -eondfi. (9) 

This mathematical problem j9]) has an unique solution^. The moment of inertia comes directly from the energy per 
unit volume of the system: E = <&tj§- + u t ■ j§[ ~ In the rotating case E = ^I ss uj 2 where I ss is the zz component 
of the inertia moment : / ss = m(n — g(n))X p f + mg(n)2 r b with T p f = J^(^^>) 2 dr, $ solution of ([9]), u, m and H 
taken to 1. This number depends on the geometry only, X r b is also a geometrical factor corresponding to rigid body 
rotational inertia (x&cy orthogonal to the axis of rotation) T r b = J n {x 2 + y 2 )dr. The relative change of the moment 
of inertia whenever the supersolid phase appears is (here I r b = mril r b) 

(I ss - I rb )/ Irb = -(1 - g{n)/n){l-l pf /l rb ) (10) 

Because T p f < T r b, one has (I ss — I r b)/Irb < as expected and observed experimentally^. The NCRI fraction 
(NCRIF) disappears, as well as the phase mode sound speed, when the supersolid phase recovers the ordinary solid 
phase (g(n) — > n). 

Within the model presented here it is easy to implement a numerical procedure to put in evidence a NCRI in a 2D 
system. We shall first minimize H — ujL z for different values of the angular frequency u, where H = ^ J | Vi/'pdr + 
\JU{t' - r)\ij}(r, t)\ 2 \ift(r\ t)\ 2 drdr' is the Hamiltonian and L z = ^ J(ip*r x V^> — ipr x Vip*)dr the angular 
momentum. The minimization should constrain a fixed total mass : N — J \ip\ 2 dr. Starting with u> = one finds the 
minimizer and then by increasing lu step by step together with the minimization procedure we follow the evolution of 
the local minima. Figure [lja represents the NCRIF as function of u>, for different values of uI/q. We observe a non- 
zero NCRIF in particular in the limit oj — > 0. Fgure[T]-6 shows this limit NCRIFo as a function of the dimensionless 
compression A = U^-jp-na? . Both curves are in qualitative agreement with recent experiments (see Fig. 3-D oftb and 
Fig. 4 o£ki). Finally, we study a gravity (or pressure) driven supersolid flow. As early suggested by Andreev et al^ 
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Fig. 1. We implement a relaxation algorithm in Fourier space with 128 x 128 modes to find a local minima in 
a square cell of 96 x 96 units, the potential range a = 4.3, for different values of Uqti. a) The NCRIF = 1 — 
L' z (uj)/ (I r b) vs. the local Maximum speed v m ax = u>L/\2 for different values of the compression parameter nllo = 
0.069, 0.084, 0.099&0.114 Here (I r b) is the average inertia moment for large ui computed numerically, b) The NCRIFo 
as a function of nUo. Note that a) and b) almost do not depend on the box size. 



an experiment of an obstacle pulled by gravity in solid helium could be a proof of supersolidity. Different versions 
of this experiment failed to show any motion^, therefore a natural question arises: How we can reconcile the NCRI 
experiment by Kim and Chan and the absence of pressure or gravity driven flows? 

In fact, our supersolid model (and it seems that supersolid helium too) reacts in different ways under a small 
external constrain such as stress, bulk force or rotation in order to satisfies the equation of motion and the boundary 
conditions. For instance, if gravity (or pressure gradient) is added then the pressure £'(n) balances the external 
"hydrostatic" pressure mgz in equation ([6]) while the elastic behavior of the solid of equation ([7]) balances the external 
uniform force per unit volume mng. No V$ nor ii are needed to satisfy the mechanical equilibria. Moreover, a flow is 
possible only if the stresses are large enough to display a plastic flow as it happens in ordinary solids. Ir£ we showed 
that a flow around an obstacle is possible only if defects are created in the crystal, in this sense we did observe a 
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plastic flow, however in the same model we observe a "supcrfluid-like" behaviour under rotation without defects in 
the crystal structure. Indeed for a small angular rotation the elastic deformations come to order uj 2 while V$ or u 
are of order oj, the equations of motion together with the boundary conditions leads to a NCRIF different from zero. 

We have realized a numerical simulation to test the possibility of a permanent gravity flow for different values of the 
dimcnsionlcss gravity Q = — ^ a . Let us consider an U-tube as in Fig[2] The system is prepared for 500 time units 
into a good quality (but not perfect) crystalline state. A vertical gravity of magnitude Q is switched-on and the system 
evolves for 500 time units more up to a new equilibrium state (see Fig. [2} a). The gravity is then tilted (with the same 
magnitude) at a given angle. A mass flow is observed at the begining from one reservoir into the other, but both 
vessels do not reach the same level eventually (see Fig. [2} b) . There is some dependence of the transferred mass on Q 
till Q w 0.0005 and the mass transfer becomes negligible from fluctuations for Q < 0.00025 indicating the existence of 
a yield-stress. The flow is allowed by dislocations and grain boundaries and it is a precursor of a microscopic plastic 
flow as in ordinary solids (e.g. ice) and as it is probably observed in Ref J^. A microscopic yield-stress could be defined 
by the smallest gravity Q such that no dislocations, defects nor grain boundaries appear. In the present model this is 
for g < 1Q- 4 . 




Fig. 2. We plot the density modulations \ip\ 2 (the dark points means a large mass concentration) of a numerical 
Simulation of eqn. |T]) with Dirichlet boundary conditions with the shape of an u-tube as in the figure. We use a 
Crank-Nicholson scheme that conserves the total energy and mass. The mesh size is dx = 1, the nonlocal interaction 
parameters are chosen as Uq = 0.01 and a = 8 (physical constants h and m arc 1), finally the initial condition is 
an uniform solution ip = 1 plus small fluctuations. One gets a crystalline state after 500 units time; then a vertical 
gravity of magnitude Q = 0.01 is switched-on, and the system evolves for 500 time unites up to a. Then gravity 
orientation is tilted in 45°. After 2000 time units the system evolves to a stationary situation b showing that the mass 
flow is only a transient. 



In conclusion, we have shown a fully explicit model of supersolid that display either solid-like behavior or superflow 
depending on the external constrain and on the boundary conditions on the reservoir wall. Our numerical simulations 
clearly show that, within the same model a nonclassical rotational inertia is observed as well a regular elastic response 
to external stress or forces without any flow of matter as in experiments^. 
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